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ABSTRACT 



Anode sheaths impact the operation of many practical plasma devices. This 
complex region is explored in detail for collisional, isothermal (identical specie 
temperatures), low-temperature plasmas, where sheath dimensions are in the micron 
range. The selected approach involves postulation of a specific electric field 
distribution with two shape factors. Previous research regarding planar anodes is 
verified and expanded upon using greater parameter ranges, 'z', a dimensionless 
quantity specifying plasma composition and condition, groups diverse plasmas into 
'families' exhibiting similar sheath characteristics, 'q', a nondimensional ratio of 
electrical energy to thermal energy in the sheatli, allows temperature effects to be 
studied. The investigation focuses on three disparate plasma families that span a z 
range of 1.1729 to 2.1493, at q values defined by plasma temperatures of 6000"K, 
3000"K, and 300‘’K. Results indicate that at lower temperatures, charge production 
in the outer sheath is generic to the electric field distribution, and that the sheaths 
themselves are nearly unaffected by substantial changes in temperature (i.e., q). 
Conversely, sheath density and extent are shown to vary significantly for differing z 
values. Newly-derived equations governing cylindrical anodes generate sheaths that 
are virtually identical to corresponding planar cases. It is shown that only those 
anodes whose radii are comparable to the plasma's 'characteristic radius' (y) must be 
treated with the cylindrical formulation; non-vacuous plasmas would require micron- 
width anodes to be thus affected. Finally, an analytical approach yields solutions that 
confirm the numerical results, and offers an algebraic approximation for high-q 
plasmas. 
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I. INTRODUCTION 



Various types of plasma devices have been in operation for over 30 years. 
Marshall's coaxial plasma gun [Ref. 1] was successfully accelerating volumes of 
hydrogen plasma as far back as the late 1950’s. In the intervening years, the number 
of applications for plasma devices has steadily increased to include 
magnetohydrodynamic (MHD) power generation, laser pumping, strategic defense, 
and electromagnetic propulsion for interplanetary spacecraft. However, due to their 
somewhat complex nature, comparatively little is understood about the sometimes- 
destructive sheath regions surrounding the electrodes in every plasma device. This 
work attempts to further previous research efforts concerning description of the 
anode sheath. 

As stated by Biblarz [Ref. 2], completely satisfactory anode sheath solutions do 
not exist for several plasma conditions: one such case involves steady, collisional, 
low-temperature, isothermal discharges. He then goes on to derive an involved, 
nonlinear differential equation that describes the entire plasma region affected by a 
planar anode, from the surface to the undisturbed plasma. A presumed (but 
'shapeable') function describing the electric field in that region, selected after much 
deliberation [Ref. 2 and Ref. 3], ultimately allows charge production rates and 
electron/ion populations to be plotted as a function of distance from the anode. 
Nondimensional parameters make the solution profiles applicable to several families 
of plasmas. 

In this work, numerical techniques are employed to verify the one previously 
solved case, and to explore planar anode solutions to several other plasma conditions. 
This can be readily accomplished due to the generality of the formulation and the ease 
with which the profiles can be produced. Sheaths and ambipolar regions are properly 
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generated for all cases. The effect of nondimensional parameter variations on the 
sheath is extensively investigated, and some general conclusions are hypothesized. In 
addition, an analytical technique is presented that supports the numerical results and 
allows for an accurate algebraic solution to low-temperature conditions. Finally, 
similar derivations produce equations that treat the cylindrical anode sheath problem; 
the equations and their profiles are then analyzed, and a comparison is made to planar 
anode findings. 
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II. THE CARTESIAN PROBLEM 



This section discusses and validates one previous anode sheath research effort 
concerning steady, low-temperature collisional plasmas. In addition, several new 
parameter cases are examined and analyzed. The geometry of the Cartesian problem 
is illustrated in Figure 1. 




Figure 1 Cartesian Anode Geometry 



A. DERIVATION REVIEW 

Previous work by Biblarz [Ref. 2] treats the planar anode sheath problem in a 
one-dimensional Cartesian fashion, and is reviewed here as it forms the starting point 
for this work. Electric field intensity (E), electron and ion population densities ( n^ 
and nj ), and all other relevant quantities are considered to vary only with linear 
distance (y) from the flat anode surface. The applicable relations are Gauss’ equation 
and the two species continuity equations, as shown on the following page. 
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(lb) 



(la) 




(Ic) 



Note that e is the electron charge constant, is the permittivity constant, the 

subscripts e and i indicate electron and ion particles, the j's are the species 
contributions to the total current density (i.e., the p's are the particle 

mobihties, and the D's are the diffusion coefficients. Although this set does not fully 
constrain the problem (as an equation describing plasma reactivity, h^, is absent), any 

selection of a specific form for E(y) implicitly fixes that variable.. It is of interest to 
explore the problem with 'guesstimates’ of E(y) in order to facilitate solutions. 

Reference 2 then introduces two new parameters, designated K* and K', which 
are defined by Equations (2a), (2b), and (2c). These K-terms are somewhat artificial 
parameters, although they are related to the total current. More meaningful is the 
derivative of K\ which is directly proportional to the net production rate of 
charges, h^ =njD- ). 




(2a) 



eD, eD 



e 



eD eD 



(2b) 



and K^==-K+ 

eD 



(2c) 
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Manipulation of the previous six equations and the Einstein relation (which 
relates the mobilities to the diffusion coefficients) yields the isothermal differential 
equation below. Note that k is the Boltzmann constant, and that all primes denote 
derivatives with respect to the variable y. 



13. 

e 






+ K^ = 






2J 

eD. 






Lv^^^oy 



I I E E^ 



(3) 



yj 



Advanced knowledge of E(y) considerably facilitates the solution of Equation (3). 
Biblarz' justification for the selected form of E(y) is detailed elsewhere [Ref. 2 and 
Ref. 3]. This function meets several critical conditions inherent to the stated 
problem, including the required distribution of E and the behavior of both K-terms at 
the boundaries. The function is reproduced as Equation (4) below. 



E(y ) = E, exp 



_(y + a)' _ 



(4) 



The parameters 'a' and 'A' are shape factors, where ’a’ is of the order of the 
sheath length, and 'A' relates to the physical parameters A(j\ and E,, (as discussed on 

pages 7 and 8). 

Substitution of the chosen E(y) into Equation (3), followed by an order-of- 
magnitude analysis, the nondimensionalization of the parameters in accordance with 
Figure 2, and further simplification based on field and current properties at the 
electrode and in the undisturbed plasma, all combine to yield the governing 
differential equation and boundary condition of Figure 3. 
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y=y/a , = (E^/Ejexp 



Uy-^iyj 



K* = 



K" 

J/^D 






V^T,y 



Z = 



^ ; 



Figure 2 Nondimensionalized Parameters 






_d_ 

(fy 



k: 

E/E, 



+ k*+ 



(mo)" 

(y + i)^ 



= 2 



K^(0) =1 



Figure 3 Simplified Nondimensional Differential Equation 
and Boundary Condition Governing the 1-D Cartesian Problem 



The introduction of the dimensionless parameters n and z further generalizes the 
problem, in that one solution to the above equation can be applicable to a large 
number of specific dimensional cases. Their additional significance is discussed in a 
subsequent section. 

As will be seen, both numerical and analytical solutions to the equations of Figure 

/ 

3 are possible, yielding profiles of and (K^) that are a function of distance from 
the planar anode. 

Further manipulation of Equations (1) and (2) derives relations that produce 
electron and ion population profiles [Ref. 2]. A more precise form of these equations 
is given in Figure 4. Note that specific solutions for are required to generate the 
n -profiles. Several distinct cases in the next section yield population curves that 
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clearly illustrate the propriety of the sheath and ambipolar regions in the plasma near 
the anode. 




Figure 4 Ion and Electron Population Equations 



Reference 2 also presents two more useful relations from analysis involving 
physical observations (Equations (5) and (6)). Recent analysis has proven that the 
infinite series in Equation (6b) converges for all values of z. These are combined and 
manipulated to yield the two important equations of Figure 5, which allow 
detennination of the parameters z and q. The procedure is as follows: 

• the type of gas defines the anode potential drop A(}\ (which is essentially 

equivalent to the ionization potential of the gas) 

• the particular case or application specifies E^ and n„ 

• z is then found by iteration of the implicit equation at the top of Figure 5 

• the choice of locks in q using the bottom equation in Figure 5 



a = 



^2e„E ^ , 

— ^ jz“ exp(2z^) 



en„ ; 

A(j), =E^^^A f(z) 



where f(z) = z 



,2m-2 



ffj(2m -l)m! 



(5) 

(6a) 

(6b) 
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Specific values of z and p allow computation of the electric field parameters A 

and a (Equations (5) and (6)), which in turn yield the tailored form of the E 

/ 

distribution (Equation 4), and ultimately permit the generation of K*, (K"^) , and , 
profiles in the plasma regions close to the anode (equations of Figure 3 and Figure 4). 

B. NUMERICALLY SOLVED CASES 
1. Procedure 

Numerical solution of the governing differential equation (Figure 3) can be 
achieved with a fourth-order Runge-Kutta FORTRAN program. First the dependent 
variable is redefined and the equations are rewritten in accordance with Figure 6. 



^[w] = 2n-q 

dy 







'(wE.r 




w- q 


L(y-^in 



w(0) = 1 



where 



w(y) = 



K-(y) 

E(y)/E, 



Figure 6 Modified 1-D Cartesian Differential Equation 
and Boundary Condition for Numerical Solution 
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Application of the Runge-Kutta scheme computes the value of w at each y . 

The data for K* is then recovered with the division of each w datapoint by the 

/ 

corresponding value for E(y)/Eo. (K*) data are extracted with the following steps 
for each value of y : 

• compute w' using the defining differential equation (Figure 6) 

• perform the operation below (derived from the Quotient Rule) 





+ K^ 



^_E^ 

Jo, 




(7) 



Finally, the h^ j profiles are computed as a function of y using the available 

information and the equations of Figure 4. A single computer routine can be made to 
perform all of the required operations (see Appendix A); the resulting data are then 
plotted. 

Note again that this technique is substantially simpler than alternate methods 
that do not presume a specific form of E(y). Ensuing discussions address the 
advantages and disadvantages of this approach. 

2. Case I: the Nitrogen Problem 

The first specific case involves verification of an earlier example [Ref. 2]. A 
planar anode in contact with nitrogen plasma was analyzed under the following 
conditions (nitrogen is a common discharge gas with well-known properties): 

• nitrogen’s anode potential drop (A4>a): 15.51 V (singly ionized) 

• the electric field strength in the undisturbed plasma (E„): 12,000. V/m 

• the particle density (n,): ICf’ m'^ 

• the temperature: 6000°K 
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The expression used to approximate f(z) (Equation 6a) was truncated after 13 
terms, producing an error in the calculation of z which is on the order of 10'*%. 

Using the given data, the following parameters result: 

• z= 1.75626 

• n = 99.12285 

• a= 1.95421x10'^ m and A = 1.17793 xlO'' m' 

• Eo = 262,265. V/m 

The computed value for E^ (obtained using Equation 4) relates specifically to 
all plasmas that are defined by the above values of z and q; this includes nitrogen 
plasma at the stated conditions. The normalized electric field E(y) for this case is 
shown in Figure 7. Note that the field decreases monotonically and abruptly to a 

constant value in the undisturbed plasma, as is required of the model. 

/ 

The corresponding K% (K*) , and h^ j profiles are presented in Figure 8 and 

Figure 9. Note that those Figures also represent aU other plasma cases whose gas 

composition, density, temperature, and electric field intensity are defined by 

z=l. 75626 and q=99. 12285. Figure 9 clearly illustrates the distinct sheath and 

ambipolar regions present for this general case.* The K-term profiles of Figure 8 

likewise exhibit expected characteristics, with K’ rising roughly monotonically from 

/ 

+ 1 to +2. while (K^) decreases to near-zero at the outer edge of the sheath. The 

/ 

leftmost downturn of (K^) , while possibly correct, is quantitatively suspect since the 

model of ionization by electron impact breaks down at the fringe of the collisionless 
region (i.e., for y/a S ~10'^). 



* The sheath is that small region, extending away from the anode, where the 
electron and ion populations are not equivalent; in the amhipolar region, the ratio 
of oppositely-charg^ particles is 1:1, but the populations of both are less than those 
that exist in the 'undisturbed plasma' 
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Electric Field (Case I) 




Figure 7 Normalized Electric Field (z=l. 75626, ri=99.1229) 
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K-term Profiles (Case I) 




Figure 8 K-term profiles (z=1.75626, q=99.l229) 
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As the parameter 'a' specific to Case I is known, it can be seen that the 
nitrogen plasma sheath for the particularly defined discharge extends approximately 
0.1 mm out from the planar anode. Similarly, the region of undisturbed plasma for 
such a device begins approximately 6 mm from the anode surface. Both of these 
results are in accordance with expectations, the generalities of which were obtained 
after much deliberation [Ref. 4]. 

These numerically-obtained Case I results confirm the work of Reference 2. 
Postulation of a suitable ionization/recombination mechanism is left to that discourse. 

The remaining objectives of this work are as follows: 

• investigation of the effects which are produced as the defining parameters are 

varied over a wide yet practical range 

• presentation of an analytical solution to the general Cartesian anode sheath 

problem 

• derivation and examination of the related Cylindrical anode sheath problem 

3. Cases II and III: Varying Temperature in the Nitrogen Problem 

Comprehensive investigation of the nitrogen plasma problem requires 
consideration of temperature's effect on the anode sheath. Thus, while all other Case 
I conditions are held constant , two somewhat more practical values for the isothermal 
temperature are considered: 

• Case II fixes TJ, at3000°K 

• Case III sets at300"K (this is representative of discharges in laser pumping) 

Temperature-only variations impact q exclusively; z and the normalized 
electric field E(y)/Eo (including the inherent factors A and a) are unaffected. Since 
q is by definition inversely proportional to (Figure 2), the decreased temperatures 
of Cases II and III signify larger values of q. As a result, the first-order derivative 
term in the governing differential equation (Figure 3) becomes less influential at the 
lower temperatures. Indeed, for the conditions of Case III, the contribution of the 
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derivative term is negligible (note, however, that this circumstance does not nullify 

the validity of the formulation; the equation does not become degenerate). 

Application of the lower equation in Figure 5 yields q values of 198.246 and 

1982.46 for Case II and Case III, respectively. The K-term profiles for these general 

cases are presented in Figures 10 and 11. Both cases exhibit the same general K-term 

tendencies that have been recognized previously, with minor variations that illustrate 

the decreasing importance of the net-ionization term. 

The effect of temperature variation (or, equivalently, q variation) on a 

plasma in the anode region is explored using Figures 12 and 13, which directly 

/ 

compare the K* and (K^) profiles of the three cases. One significant observation. 

contributed by Biblarz subsequent to Reference 2, hypothesizes that the asymptotic 
behavior of the K-terms with decreasing temperature reflects the approach to 
ionization that is generic to the electric field distribution; the overall results are thus 
independent of the details of ionization and recombination. As an important practical 
example, measurements of charge density would not reflect relevant details of the 
ionization/recombination mechanisms. 

Conversely, the h^; profiles for this general z=l. 75626 case are nearly 
unaffected by the changes in temperature. Figure 14 illustrates charged particle 
curves that exhibit only a barely perceptible shift of position over the entire range of 
temperatures tested. This unexpected result has substantial implications; the most 
significant of these is the acceptability of the isothermal-particle assumption (i.e., the 
assertion that % =Tj =To). Though such a premise may be unrealistic, its validity 

appears to be inconsequential to the population profiles. Electric-field effects 
dominate any temperature effects in the cases examined. In fact, the physical 
significance of q (stated in Reference 2 as a ratio of electrical to thermal energy) 
predicts such dominance with its large magnitudes. 
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K-term Profiles (Case II) 




Figure 10 K-terni profiles (z=1.75626, r)=198.246) 
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K-term Profiles (Case III) 




Figure 11 K-term profiles (z=1.75626, q=1982.46) 
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Kplus Profiles for Varying Temperatures 




Figure 12 K* profiles (z=l. 75626) 
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(Kplus)^ Profiles for Varying Temperature 




Figure 13 (K^)' profiles (z=1.75626) 
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n-Profiles for Varying Temperature 




Figure 14 n^., profiles (z=l. 75626) 
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C. EXPLORING THE ENVELOPE FOR z AND n 

The foregoing results include a preliminary investigation into the effect of r)- 
variation on the planar anode sheath. This upcoming section more fully explores how 
changes in the parameters z and x] impact the sheath region. 

1. Practical Extremes of z 

The parameter z is a dimensionless representation of three basic variables 
(A(]), , E» , and n„) that specify the condition of a plasma; several different gases, at 

differing densities and field strengths, can be related by their equivalent z values. 
Various diverse plasmas belonging to the same 'z-family', at the same temperature, 
would theoretically exhibit identical sheath characteristics. Using the technique 
discussed previously, anode sheaths for plasmas defined by two limiting values of z 
are examined. 

Singly-ionized elements possess anode potential drops ( A(J)^ ) that range from 

24.46 Volts (helium) down to 3.87 Volts (cesium), with the low values being 
preferred. Common densities (n^) for these type collisional plasmas var)' 

approximately from ICy* to 10‘° particles per cubic meter. Finally, reasonable field 
strengths in the undisturbed plasma (E„) are somewhat arbitrarily set from 5000 to 

30,000 Volts per meter. Application of the upper equation in Figure 5 yields the 
following 'practical' extremes for z: 

1.1729 s z s 2.1493 (8) 

This relatively small numerical range represents an extremely diverse 
grouping of plasmas; z is a decidedly sensitive parameter. Small z values represent 
those plasmas with small potential drops and low densities that are subjected to very 
high electric fields (e.g., low density cesium at 30,000. V/m). Plasmas with large z 
values are characterized by gases with large potential drops, at conditions of high 
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density and low field strength (e.g., high pressure helium at 5000. V/m). By way of 
comparison, the nitrogen plasma considered earlier in this work (median potential 
drop and density, small field strength) has a slightly-larger-than-median z value of 
1.7563. 

Values for r| are required to completely constrain the problem, which entails 
the choosing of one or more temperatures. To facilitate direct comparisons, the same 
arbitrary quantities for that have been considered for the nitrogen case are 

designated as the standard values (i.e., 6000°K, 3000°K, and 300°K). 

Consequently, numerical sheath solutions are next generated and examined 

for these forementioned hmiting conditions: 

• 'Small z' plasmas (z= 1.1729) at the three standard temperatures (q= 16.4135, 

32.8269, and 328.269 for T; at 6(X)0°K, 3(XK)°K, and 300°K) 

• 'Large z' plasmas (z= 2.1493) at the standard temperatures (q= 257.563, 

515.125, and 5151.25 as T; decreases to300°K) 

2. Sheath Solutions for 'Small z' Plasmas 

K-term profiles for the 'small z' plasmas (at the three standard temperatures) 
are depicted in Figures 15, 16, and 17. Qualitatively, each of these curves displays 
the same previously noted characteristics inherent in the nitrogen plasma curves: 
monotonically increasing and well-behaved K*, and a rate-of-charge-production term 

(K^) that decreases and vanishes at the outer edge of the sheath. As before, the 

/ 

leftmost downturn of (K*) is quantitatively suspect at the fringe of the collisionless 
region. 

These reassuring results extend to the composite curves of Figures 18 and 19, 
which depict K-term profile variations as a function of temperature (i.e., q) changes. 
In particular, the charge production rate curves of Figure 19 again illustrate the 
striking imphcation that charge production in the outer sheath is independent of the 
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K-term Profiles, Low z, T=6000°K 




Figure 15 K-term profiles (z=1.1729, r|=16.4135) 
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K-term Profiles, Low z, T=3000'’K 




Figure 16 K-term profiles (z=1.1729, r|=32.8269) 
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K-term Profiles, Low z, T=300°K 




Figure 17 K-term profiles (z=1.1729, r)=328.269) 
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Kplus Profiles, Low z, Various Temp. 




Figure 18 K* profiles (z=1.1729) 
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(Kplus) Profiles, Low z, Various Temp. 




Figure 19 (K^) profiles (z=1.1729) 
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temperature. Moreover, the inner sheath once again sees an increase in the net 
production of charges as the temperature decreases. 

Electron and ion population profiles for the 'small z' plasmas are given in 
Figure 20. As with the nitrogen case (and other plasmas of that z-family), the sheath 
and ambipolar regions are clearly visible and quite appropriate, while the population 
curves themselves again appear to be nearly unaffected by large variations in 
temperature. 

Such qualitative observations for the 'small z' plasmas indicate that previous 
impressions concerning 'median z' plasmas (e.g., nitrogen) are not limited to that one 
category, but are possibly characteristic of ah collisional low-temperature plasmas. 
The sheath profiles generated for the 'large z' cases of the next section further 
confirm this hypothesis. 

Quantitative deviations in the sheath due to the change of z are addressed in a 
later section. 

3. Sheath Solutions for 'Large z' Plasmas 

The high density, small field 'large z' plasmas continue the trends established 
by the other collisional plasmas, albeit with some minor differences. One such 
difference is illustrated by the K-term profiles of Figures 21, 22, and 23. The curves 
display the same general traits of the other plasmas, with the exception that the 
curve no longer climbs monotonically. This small anomaly becomes more 
pronounced at higher temperatures (Figure 21), yet does not appear to significantly 
impact the charge production rate curves, which still behave as expected. 

Similarly, there are no surprises in the composite curves of Figures 24 and 
/ 

25. In fact, the (K*) curves of Figure 25 demonstrate even more strikingly the 
generic, electric-field-dependent nature of ionization and recombination in the sheath. 
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n-Profiles, Low z. Varying Temp. 




Figure 20 ifi^, profiles (z=1.1729) 
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K-tei-m Profiles, High z, T=6000°K 




Figure 21 K-term profiles (z=2.1493, ri=257.563) 
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K-term Profiles, High z, T=3000°K 




Figure 22 K-term profiles (z=2.1493, r|=515.125) 
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K-term Profiles, High z, T=300°K 




Figure 23 K-term profiles (z=2.1493, r|=5151.25) 
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Kplus Profiles, High z, Various Temp. 




Figure 24 K* profiles (z=2.1493) 
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(Kplus/ Profiles, High z, Various Temp. 




Figure 25 (K*) profiles (z=2.1493) 



The n-profiles of this class of plasmas are depicted in Figure 26. While there 
are some yet-to-be-addressed differences from previous n-profiles, the sheath and 
ambipolar regions are again evident, as are the profiles' theorized independence from 
temperature changes. 

4. The Effect of z on the Planar Anode Sheath 

Results to this point have revealed some important and perhaps non-intuitive 
concepts concerning the properties of planar anode sheaths, concepts which appear 
applicable to aU low-temperature, collisional plasmas. However, sheath variations 
attributable solely to changes in the parameter z have yet to be discussed. 

Compeuison of the n-profiles for each z-family (Figures 14, 20, and 26) 
illustrates a subtle point: although all n-profiles have been shown to be nearly 
independent of temperature, they do become slightly more sensitive to temperature 
variations as the value of z decreases. This is possibly attributable, in part, to the 
lower densities of 'small z' plasmas; the temperature- dependent kinetic energy 
changes of the particles may produce a more observable effect on the total 
ionization/recombination process in a plasma not already saturated with density- 
driven collisions. Perhaps more contributory is the previously-discussed dominance 
of electrical energy over thermal effects. The 'small z' plasmas exhibit smaller 
magnitudes of than other plasma families, decreasing such dominance. The 
involvement of Act)^ in this phenomenon remains unclear. In any event, such n- 
profile temperature shifts remain nearly negligible for all plasmas considered to this 
point, which further validates the innocuous nature of the constant temperature 
assumption. 

Figure 27 depicts the curves of three diverse z-families at the same 
temperature (6000° K). As noted previously, all are well-behaved and do not differ to 
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n-Profiles, High z, Varying Temp. 




Figure 26 n^ j profiles (z=2.1493) 
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Kplus Profiles, 6000°K, Various z 




Figure 27 at6000’’K for Different z-Families 
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any significant degree. Contrastingly, the charge production rate curves of these 
three plasma families (Figure 28) are somewhat disparate and very insightful. 'Large 
z' plasmas show significant charge production throughout the entire sheath, in 
contrast to the localized and decreased production rates depicted for the 'small z' 
plasmas. This may indicate the fact that electric-field effects are more localized for 
smaller z. 

The n-profiles for the three z-families are depicted in Figure 29. Unlike 
similar profiles for varying temperature (i.e., p), these curves are decidedly affected 
by changes in z. Most prominendy, the magnitude of the charged particles changes 
dramatically as z varies, especially in the sheath. This is a graphic indication of z- 
induced changes in both Eq and current flow for diverse plasmas. Also noteworthy is 

the fact that the sheath itself extends further from the anode surface as the value of z 
decreases. In particular, the sheath of the 'small z' plasmas stretches over ten times 
the distance occupied by the 'large z' sheath. 

5. Some Final Thoughts on the Influence of r| 

The effect of r| variations on the anode sheath has already been indirectly 
addressed via the extensive consideration of temperature's influence on the problem. 
The bottom equation of Figure 5 inversely relates these two parameters for fixed z 
and A(})^ . Physical interpretation of q as a ratio of electrical and thermal energy has 

also been reviewed. 

The purely numerical importance of q is apparent with a glance at the 
governing differential equation (Figure 3); its reciprocal controls the influence of the 
charge production rate term. Large values of q reduce this term to near zero, the 
consequences of which can be seen in every low-temperature case (i.e., 300"K). To 
fully and completely explore this parameter's effect on the sheath, solutions are 
presented for the original nitrogen plasma at artificially small values of q (in 
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(Kplus)-" Profiles, 6000°K, Various z 




Figure 28 (K*)* at 6000"K for Different z-Families 
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n-Profiles as a Function of z, T=600CfK 




Figure 29 n,,, Profiles at 6000"K for Different z-Families 
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particular, q is set equal to 25., 10., and 1.0). Note, however, that the validity of 

these curves is questionable at best, since such small q values equate to some 

extraordinarily high temperatures for the nitrogen (up to 595,000”K!), the existence 

of which contradicts the original assumption of low-temperature plasmas. 

Nevertheless, the following curves may peiliaps reveal some valid trends. 

Despite the unusual and probably erroneous K* plots for the small-q 

/ 

conditions (Figure 30), the (K*) curves of Figure 31 hearteningly retain some of the 
familiar traits concerning charge production rate in the sheath. The corresponding n- 
pro files are presented in Figure 32; remarkably, only the anomalous 'q=l' profiles 
exhibit tendencies somewhat contrary to those already observed. 

D. ANALYTICAL SOLUTIONS 

In order to confirm the validity of the numerically-obtained solutions, an 
analytical solution to the problem's governing relation (Figure 6) is presented. 

1. Procedure for 'The Outer Expansion Method' 

Reference 5 offers an approach which is ideally suited to the form of the 
differential equation in Figure 6. The tenn (l/q) is justly defined as a 'small' 

parameter, and the equation is already in dimensionless form. A series solution for w 
is presumed, approximate yet valid for 'non-small' values of y . Equation (9) depicts 
the power series expansion for w: 



w 



-I 



W, 



k=0 



'"-T 

n. 



= Wo + 



'-1 



Wi + 



fiY 



w, + 






Wj +••• 



(9) 



Substitution of Equation (9) into the governing equation allows terms with 
like powers of (l/q) to be equated, which in turn allows each w,; to be solved for 
analytically. Because the solution for w,, is purely algebraic in this formulation, the 
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Kplus Profiles, z for Nitrogen, Low eta’s 




Figure 30 K* profiles (z=1.75626, Small n'*) 
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(Kplus)^ Profiles, z for Nitrogen, Low eta’s 




Figure 31 (K')' profiles (z=1.75626, Small ti*^) 
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n-Profiles, z for Nitrogen, Low eta’s 




Figure 32 j profiles (z=1.75626, Small n*^) 
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first-order derivatives in each of the subsequent expressions can be evaluated 
exactly (i.e., symbolically). The total solution w is thus an infinite summation of 
exact expressions. Prevalent magnitudes of q allow highly precise approximations of 
w with the series truncated to only three terms. 

Utilization of this technique produces the following approximate analytical 
'outside' solution (wj to the governing differential equation of Figure 6: 



W' = W = Wn + 
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Figure 33 Three-Term Approximate Analytical Solution (Wa) 



The structure of the equations themselves offers insight. Note that each w,^ is 
a function only of z and y ; they are completely independent of q (i.e., temperature). 

The contribution of q to the total solution manifests itself solely in the series- 
expanded expression for w^. Note also that the equation for w,, is identical to the 

original governing differential equation for the case when q is allowed to increase to 
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infinity. As a result, the partial solution Wq is itself a valid approximation for the 

true solution (w) at large values of r|. This fact is also corroborated by the series- 
expanded expression for w^; all terms beside vanish for large r|. 

As before, plots of K’ can be recovered from the data for w (or, in this case, 
w^). Using parameter values from the first nitrogen plasma case (z= 1.7563 and 

n=99.1229), profiles from both the 'Outer Solution' and the Fourth-Order Runge- 
Kutta scheme are compared (Figure 34). Although the 'Outer Solution' does indeed 
appear to diverge for extremely small values of y, the two curves are nearly 
identical. Such correlation offers comforting evidence for the validity of the 
numerical procedures and results presented in this work. 

Figure 35 illustrates the soundness of using just w^ to recover data for large- 

r) cases. The solid curve is from the Runge-Kutta solution for the 300“K nitrogen 
case (z=1.7563 and n=l 982.5); the dotted profile depicts data generated using only 
the partial solution w^. 
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Kplus (Case I): Numerical vs. Analytical 




Figure 34 Comparison of Numerical and Analytical Solutions 
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Kplus (High eta): Numerical vs. Analytical 




Figure 35 One Term (u„) Algebraic Approximation for High q's 
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III. THE CYLINDRICAL PROBLEM 



The following sections present detailed derivations and analysis for the 
cylindrical anode sheath problem. As with the planar cases, the focus is limited to 
steady, low-temperature collisional plasmas. Exploration of the cylindrical sheath is 
desirable, as several functioning plasma devices contain or employ a cylindrical 
architecture. In addition, these results can be compared and contrasted to the 
previously-presented Cartesian results. The geometry of the one-dimensional 
cylindrical problem is illustrated in Figure 36. 



Attainment of a governing differential equation for the cyUndrical case begins 
with the re-derivation of the applicable relations (Gauss' equation and the two species 
continuity equations) in cylindrical form. Reference 6 presents the following 
differential form of Gauss' law; 




(distance from 
the center) 



r 



Figure 36 Cylindrical Anode Geometry 



A. DERIVATIONS 




( 10 ) 
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V is the del operator, E is the vector form of the electric field, p is the charge 
density (which can be represented here in accordance with Equation 11), and £„ is the 

permittivity of free space. The divergence (V E) of the electric field in cylindrical 
coordinates is given in Equation (12) below: 

p = e[n, -nj (11) 



V E = 




dE 

+ 2 - 

dz 



( 12 ) 



For a one-dimensional electric field that varies only as a function of radius (r), 
combination of Equations (10), (11), and (12) yields the following one-dimensional 
cylindrical form of Gauss' equation : 



dr 




(13) 



Comparison of Equations (13) and (la) accentuates the appearance of a new term 
which is the result of the cylindrical derivation. 

The species continuity equations are derived from conservation equations [Ref. 7] 
that can be manipulated into the form given below: 

-j, = -ep;i,E-e4Vn^ (14a) 

ji=ePiniE-eD,Vnj (14b) 



For a one-dimensional electric field, the previous relations are easily rewritten as 
shown in Equations (15a) and (15b) on the following page; these are the 1-D 
cylindrical species continuity equations . 
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j.=eM.n,E + eD^^ 
dr 

dn 

J. =eM.niE-eD,— ^ 



(15a) 

(15b) 



Note that these species equations are nearly identical to those for tlie planar case 
(Equations lb and Ic); only the independent variable has been renamed. 

The procedure from this point is the same as that of Reference 2: 

• Combine these forms of Gauss’ equation and the species continuity equations 

with the Einstein relation 

• Incorporate the K-teims as defined by Equation (2) 

• Assume isothermal conditions 

• Derive a single nonlinear differential equation in K\ E, and their derivatives 
The desired intermediate differential equation (the cylindrical counterpart to 

Equation 3) is presented below. Note that all primes here denote derivatives with 
respect to the variable r. 
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(16) 



In compjirison to Equation (3), several new terms exist that are due solely to the 
cylindrical nature of the derivation. However, the validity of the above expression is 
demonstrated by letting r increase towards infinity, which causes the equation to 
approach the planar case. In that limit. Equation (16) reduces exactly to the 
corresponding Cartesian expression of Equation (3). 

The next step calls for an 'educated guess’ of the form for E(r), and substitution 
of that function into Equation (16). For reasons stated previously, the basic form of 
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the electric-field function given in Equation (4) will be retained, although in a form 
modified for use with a cylindrical anode. The plasma electric field begins at the 
anode surface, which for a cylindrical anode is at r=i|) (as opposed to r=0, the anode 

center). Thus, the presumed electric field function is 



E(r) = E, exp 



B 



L((r-r.)+b)-J 



(17) 



where i;, is the anode radius, B and b are constants that give specific shape to the 
field, and the range of r is from i;, to °°. 

Following substitution of E(r) into Equation 16, the resulting lengthy expression 
is made dimensionless through use of the relations and definitions given in Figure 37. 
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Figure 37 Cylindrical Nondimensionalized Parameters 



Note that q is the cylindrical counterpart to z, and that r has a more conventional 
range of zero to infinity. The new parameter y appears because the simplifications 
that were previously employed (based on field and current properties at the electrode 
surface and in the undisturbed plasma) in the Cartesian derivation [Ref. 2] cannot 
eliminate all of the same constants when applied to the cylindrical case, y has units of 
length, and can be thought of as a 'characteristic radius' of the plasma conditions and 
anode width that define q and q^. 
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Further simplification of the differential equation is achieved with the derived 
expression [Ref. 2]: 



J ] 




fen,E^] 






1 j 



(18) 



Finally, order-of-magnitude analysis reveals that some of the surviving terms 
cannot be neglected unless the anode radius (i;,) is sufficiently large. As a result, the 

differential equation can take one of two possible forms; the appropriateness of either 
is a function mainly of the anode’s radial magnitude. Thus, the governing differential 
equation and boundary condition for 'normal-sized' anodes (e.g., mm) are 

given in Figure 38, while the corresponding equations for a 'wire-thin' anode (e.g., 
i;, «0.1 mm) are given in Figure 39. 




Figure 38 Simplified Nondimensional Differential Equation 
and Boundary Condition, Governing the 1-D 'Normal Anode' 

Cylindrical Problem 
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Figure 39 Simplified Nondimensional Differential Equation 
and Boundary Condition, Governing the 1-D 'Wire Anode' 

Cylindrical Problem 
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The equations of Figure 39 are actually appropriate for aU magnitudes of i;,, but 
some of the terms become negligible for the larger anodes. In fact, the degree of 
influence exerted by the extra terms of that differential equation (Figure 39) is 
directly tied to the ratio y/r^ ; unless the cylindrical anode’s radius is smaller than, or 

within an order-of-magnitude of, the 'characteristic radius’ (y), then those extra 
terms are insignificant For this reason, the simplified differential equation of Figure 
38 is also provided. 

Note also that this ’normal anode’ differential equation is identical in form to the 
corresponding planar equation (Figure 3), which indicates that in spite of a tortuous 
derivation, most cylindrical and planar anodes disturb the plasma in nearly the same 
manner. Finally, as ij, in the expression of Figure 39 increases toward infinity 

(approaching the planar case), the extra terms vanish and the planar differential 
equation is recovered. 

The parameters q and q,, have comparable meaning to their planar counterparts, 
and their specific values are computed using similar techniques. It can be shown 
through extensive manipulations that the cylindrical equivalents to Equations (5) and 
(6) are as follows: 




(19) 



A4), =E^^f(q) 



(20a) 




(20b) 
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Note that the infinite series converges for all values of q, and that the y/fo 
in tlie denominator of Equation (19) is only relevant for the 'wire-thin' (and smaller) 
anodes. The above equations are combined to generate the two important equations 
of Figure 40, which are used to compute the parameters q and in accordance with 



the following familiar procedure: 

• the gas composition defines 

• the particular case or application specifies E„ and n^ 

• the device in use fixes the anode radius ij, 

• q is then computed using the implicit equation at the top of Figure 40 

• the plasma temperature TJ, then yields via the bottom equation of Figure 40 

• for the extremely thin anodes, y is computed Irom its definition and the value of 

E(f = 0); both expressions are found in Figure 37. 
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Figure 40 Equations That Yield q and q^ 



As with the planar cases, specific values of q and q^. allow computation of the 

electric field parameters B and b (Equations 19 and 20), which in turn yield the fitted 

form of the E distribution (Equation 17), and ultimately permit the creation of K*, 
/ 

(K^) , and h^ i profiles in the plasma regions close to the anode. All that remains to 
be derived are the cylindrical forms of the equations that compute those n-profiles, 
which are produced through manipulation of Equations (2), (13), and (15). The 
resulting elaborate expressions are presented in Figures 41 and 42. 
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(r + 1) (r + ^,/b)J 

Figure 42 Electron Population Equation 

Gratifyingly, these n-profile equations produce their planar counterparts as i;, 
approaches infinity. Note also that the expressions for fij and n^ are valid for all 
cylindrical anodes. Unfortunately, only one term in Figure 41 could be neglected for 
larger, 'normal-sized' anodes, which stifles the notion of a separate 'simplified' 
equation for such cases. 

B. NUMERICALLY SOLVED CASES 
1. Procedure 

The previous derivations allow detailed numerical analysis of cylindrical 
anode sheaths in certain plasma conditions, using a procedure similar to that 
performed on the planar problem. As before, the governing equation is rewritten in 
terms of w, as defined in Equation (21). The resulting differential equations and 
boundary conditions for both the 'normal-sized' anode (~5 mm radius and larger; this 
form also applies when [yf r,)« 1) and the 'wire-thin' anode (-0.5 mm radius and 
smaller, or when (Y/ro)“0 cases are presented in Figures 43 and 44. Note that the 
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form of the equations for the 'normal-sized' anode is identical to that of the 
corresponding planar equations (Figure 6). 



w( f ) = 
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E('r)/E, 



( 21 ) 
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Figure 43 Modified 1-D Cylindrical Equations 
for Numerical Solution ('Normal-Sized' Anode) 
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Figure 44 Modified 1-D Cylindrical Equations 
for Numerical Solution ('Wire-Thin' Anode) 



As before, application of a FORTRAN Runge-Kutta algorithm (see Appendix 

A), combined with the previously-outlined data manipulations, yield the numerical 
/ 

profiles for K\ (K^) , and j for any desired plasma and anode width. 

2. The Nitrogen Problem with a Cylindrical Anode 

The nitrogen plasma conditions designated previously as 'Case I' for the 
planar anode are now analyzed for both the 'normal-sized' and 'wire-thin' cyhndrical 
anodes. The Case I conditions are repeated on the next page: 
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• nitrogen’s anode potential drop (A4>,); 15.51 V (singly ionized) 

• the electric field strength in the undisturbed plasma (E, ); 12,000. V/m 

• the particle density (n„): lO’’’ m‘^ 

• the temperature (TJ,): 6000“K 

• the anode radii under consideration (ij)): 10.0 mm and 0.1mm 

These specifications and the newly-derived equations produce the following 
parameter values for their respective anodes: 



TABLE 1 CYLINDRICAL PARAMETERS FOR CASE I CONDITIONS 




Y 



b - 



1 ^ = 10. mm 
1.75656 
99.1571 
N/A 

1.953M0'^ m 
262,543. V/m 



ij, =0.1. mm 
1.78774 
102.818 
3.9597 xlO - m 
1.813x10'^ m 
293,216. V/m 



A comparison reveals that the parameter values for the 'normal-sized' (10. 
mm) anode differ negligibly (only 0.1% and less) from their corresponding planar 
values. This is not an unexpected result, in light of the cylindrical equations' 
similarity to their planar counterparts for the bigger anodes. More surprising are the 

values computed for the 'wire-thin' (0. 1 mm) anode; while their differences from the 
planar values are not insignificant, only differs by more than 8%. The extreme 

sensitivity of the parameter q, however, implies that percentage comparisons may be 
misleading. 
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Direct comparison of tlie planar and cylindrical profiles best illustrates the 
near-negligible differences in their respective sheath and ambipolar regions. Figures 
45 and 46 present the Case I K-terin profiles produced by the three anode types 
(plane, 'normal' (10. mm) cylinder, and 'wire' (0.1 mm) cylinder), while Figure 47 
depicts the n-profiles for each of those anodes. As can be seen, profiles for the 
'normal' cylindrical anode are visually indistinguishable from the planar anode 
profiles. Amazingly, even the 'wire' anode profiles are nearly identical to those of 
the planar anode. This is a welcome result; despite many differences in the derivation 
and appearance of the governing equations for each type of anode, it seems that for 
this 'small y' family of plasmas those cylindrical anodes whose radius is greater than 
that of a human hair can be treated to behave as planar . For such plasmas, only those 
cylindrical anodes whose radii are on the order of the sheath thickness itself would 
exhibit significant non-planar characteristics. 

To further validate this hypothesis, the three anodes are also compared under 

the lower-temperature 'Case III' circumstances, which represent the same conditions 

as Case I except that = 300"K. Consequently, the only parameter listed in Table 1 

that changes is q,.; its Case III values are 1983.14 for the 'normal-sized' cylinder, and 

/ 

2056.36 for the 'wire' anode. Figures 48, 49, and 50 show the K% (K^) , and j 

profiles (respectively) for Case HI conditions. Once again, the planar and 10-mm- 
radius anodes disturb the plasma in identical fashion, while the 0. 1 -mm-radius anode 
exhibits nearly negligible variations. 

It is worth noting that the ratio (v/ro) can theoretically become a meaningful 
factor even for the 'normal-sized' anodes; extremely low-density plasmas could 
possibly generate values of y in the millimeter range. Thus, the notion of sheath 
characteristics being strictly a function of anode radius is valid only for the 
conditions being considered in this work. 
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Kplus, Cyl. Anode Comparison, Case I 




Figure 45 Profiles, a Three-Anode Comparison 
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(Kplus/, Cyl. Anode Comparison, Case I 




Figure 46 (K*)’ Profiles, a Three-Anode Comparison 
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n-Profiles, Cyl. Anode Comparison, Case I 




Figure 47 n^, Profiles, a Three-Anode Comparison 
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Kplus, Cyl. Anode Comparison, Case III 




Figure 48 K* Profiles, a Three-Anode Comparison 
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(Kplus)^, Cyl. Anode Comparison, Case III 




Figure 49 (K*)' Profiles, a Three-Anode Comparison 
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n-Profiles, Cyl. Anode Comparison, Case III 




Figure 50 n,., Profiles, a Three-Anode Comparison 
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IV, CONCLUSIONS 



The manner in which one-dimensional planar and cylindrical anodes disturb 
collisional, isothermal plasmas has been explored in detail. Appropriate sheath and 
ambipolar regions have been shown to develop, for a given specific and reasonable 
electric field distribution, across a wide range of plasma types and conditions. 
Numerical solutions of the nondimensional governing equations have allowed specific 

and important observations to be made concerning; 

• the effect of temperature variation on plasmas in the anode region 

• traits that are characteristic of various plasma 'families', and traits that appear 

to be common to all collisional, low-temperature plasmas 

• the effect that differing anode radii have on the sheath for cylindrical anodes 

• the degree of difference betw'een sheaths of planar and cylindrical anodes 

In addition, an analytical method has been presented that offers a simplified yet 

accurate algebraic solution for the lower-temperature plasmas. . 

The following briefly summarizes the specifics for each of these results. 

/ 

A ppropriate Sheath Existence : In each case considered, the (K*) curves showed 

that the net rate-of-charge-production decreases and vanishes toward the outer edge 
of the sheath. Moreover, each of the n-profiles clearly depicted sheath and ambipolar 
regions which were reasonably located with respect to the anode surface. 

Overall-Temperature Variation Effects : Results of this work indicate that at 
lower temperatures, charge production in the outer sheath is generic to the electric 
field distribution, and is thus independent of the details of ionization and 
recombination. Additionally, the ; profiles, and thus the sheaths themselves, have 

been shown to be nearly unaffected by substantial changes in temperature. This 
somewhat unexpected result, which held for every plasma case considered, reduces 
any consequences caused by the isothermal-particle assumption (i.e., the assertion that 
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= T; = Ty). It is interesting to note tliat the large magnitudes and physical meaning 
of r) (the ratio of electrical energy to thermal energy in the sheath) actually predict 
such independence from temperature variations. 

Plasma 'Family' Characteristics : All lo\v-temperature, collisional plasmas can be 
partitioned into broad groupings based on their specific nondimensional 'z' or 'q' 
value (z denoting a planar anode problem, and q signifying cylindrical anodes). Each 
plasma 'family' thus exhibits its own characteristic electric field and sheath 
properties, even though its members may be widely diverse in composition and 
condition. 

Plasmas with large z (or q) values show significant charge production throughout 
the entire sheath, in contrast to the localized and decreased production rates depicted 
for the 'small z/q' plasmas. This may indicate the fact that electric-field effects are 
more modest and thus localized in those plasmas with smaller values of z or q. In 
addition, the magnitude of the charged particles near the anode changes dramatically 
as 'z/q' varies, especially in the sheath. The sheaths of 'large z/q' plasmas exhibit 
charged-particle densities up to one order-of-magnitude smaller than those of the 
'small z/q' sheaths. More importantly, the sheath itself extends further from the 
anode surface as the value of 'z/q' decreases. In particular, the sheath of the 'small 
z/q' plasmas stretches over ten times the distance occupied by the 'large z/q' sheath. 

Cylindrical Anode Radius Effects : The sheath and the charge production rates 
have been shown to vary as the anode radius varies; however, the differences are 
almost negligible for all practical anodes and conditions. It appears that only 
extremely low-density plasmas, or cylindrical anodes whose radii are on the order of 
the sheath thickness itself, would exhibit significantly different characteristics. 
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Planar and Cylindrical Anode Differences : In an extension of the previous 
paragraph, the results of this work indicate that for the range of conditions explored, 
all practical cylindrical anodes can be treated to behave as planar anodes. 

Simple Analytical Approximations (Low Temperatures Only) : Using the 'Outer 
Solution’ technique, an algebraic solution can be used to generate highly accurate 
approximations to low-temperature anode problems. In addition, as long as 
derivatives of the presumed electric field distribution exist, the sheath profiles for all 
cases can be produced analytically to any desired accuracy. 

One suggestion for further work involves researching the effects that other 
electric field distributions would have on the anode problem. Such distributions may 
be derived from empirical data or from other 'guesstimates', and may have three or 
more adjustable parameters. The most prominent reason for this suggestion is the 
need to corroborate the generic charge production rate observed at lower 
temperatures. 

Another area requiring further study is the exploration of two-temperature 
plasmas and their anode sheaths. Results for such conditions could confirm the 
assertion of sheath insensitivity to temperature changes. 
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APPENDIX A 



Applicable FORTRAN Programs 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



DOUBLE PRECISION YA ( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 
DOUBLE PRECISION XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 
REAL PI ,XL,XCHK, EINF, EO , Z , ETA 
INTEGER NS, PI 
PRINT * 

PRINT FOURTH ORDER RUNGE-KUTTA SCHEME ' 

PRINT FOR THE CARTESIAN COMPUTATION ' 

PRINT OF W, KPLUS, & KPLSPRM' 

PRINT *, 'INPUT VALUES FOR EINF, EO, Z, AND ETA' 

READ *, EINF, E0,Z, ETA 

IM=1 ! Number of equations 

Y(l) = 1.00 ! Initial condition for y~l at x=XP . 

Y(2) = 0 ! Initial condition for y~2 at x=XP (if nec) 

PRINT * 

PRINT *, 'INTERVAL OF X FOR PRINTING ?' 

READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-x)' 

READ *,H 

NS = NINT(PI/H) 

PRINT *, 'MAXIMUM X TO STOP CALCULATION ?' 

READ *,XL 

PRINT * , ' H= ' ,H 

Pl = 0 
XP=0 
HH=H/2 
KPLUS=1 . 0 



print *, ' HH= ' ,HH 
print * , ' NS= ' , NS 
PRINT * 



LI = 0 ! Line no. initialization 

PRINT*, 'LINE y/a KPLUS' 

WRITE (*,98) LI,XP, KPLUS 



LI=LI+1 
DO N=1,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 



! Old time 
! New time 
! Midpoint time 



J=1 ! This part computes k~l. 

DO 1=1, IM 

YA(I)=Y(I) 

END DO 
XA=XB 
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on no oo oo 



CALL FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA ) 

J=2 ! This part computes k~2. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END DO 
XA=XM 

CALL FUNCT(EK, J,YA,H,XA, EINF, E0,Z, ETA) 

J=3 ! This part computes k~3. 

DO 1 = 1, IM 

YA( I )=Y( I )+EK( 2 , I )/2 
END DO 
XA=XM 

CALL FUNCT( EK, J, YA,H,XA, EINF , EO , Z , ETA ) 

J=4 ! This part computes k~4. 

DO 1 = 1, IM 

YA( I )=Y( I )+EK( 3, I ) 

END DO 
XA=XP 

CALL FUNCT(EK, J,YA,H,XA, EINF, E0,Z, ETA) 

DO 1=1, IM ! 4-th order Runge-Kutta scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2, I ) *2+EK( 3, I ) *2+EK( 4 , I ) )/6 
END DO 



c 



W=Y( 1 ) 

You now have W( y/a ) -- to get K+( y/a ) you must 
multiply by { E( y/a ) / E(0) } 

E=(EINF/EO)*DEXP( (Z**2. )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

To get the derivative of W [ W' ( y/a ) ] , just plug the 

computed values of W back in the original Ist-order ODE 

J=5 

YA( 1 )=W 
XA=XP 

CALL FUNCT( EK,J,YA,H, XA, EINF, E0,Z, ETA) 

WPRM=EK(5,1)/H 

NOW, to get d/dy [ K+ ] from d/dy [ K+/(E/E0) ] , must 

perform the following operation: 

KPLSPRM=(WPRM*(E**2. ) +KPLUS*E* ( -2 . * ( Z**2 . )/( (XP+1. )**3. ) ) )/E 

To keep from generating unplottable 50,000 point data files, 
the following will edit out data points depending on 
'where' they occur 

XCHK=(XP/H) 

IF (XCHK .GT. 100000) GOTO 72 
IF (XCHK .GT. 10000) GOTO 73 
IF (XCHK .GT. 1000) GOTO 74 
IF (XCHK .GT. 100) GOTO 75 
GOTO 44 
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Pl=Pl+l 
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10000) GOTO 88 



IF (Pi .NE. 

Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi .NE. 1000) GOTO 88 
Pl = 0 
GOTO 44 
Pl=Pl+l 

IF (Pi .NE. 100) GOTO 88 
Pl = 0 
GOTO 44 
Pl=Pl+l 

IF (Pi .NE. 10) GOTO 88 
Pl = 0 



WRITE (13,*) XP, KPLSPRM 
WRITE (12,*) XP, KPLUS 
WRITE (10,*) XP, W 

CONTINUE 



END DO 

WRITE (*,98) LI,XP, KPLUS 
FORMAT(lX, I2, FIO.6, 2X, 1 p4e16.8) 

IF (XP .LT. XL) GOTO 28 

PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ *,K 

IF(K.EQ.l) GOTO 1 
PRINT* 

END 






SUBROUTINE FUNCT ( EK , J , YA , H , XA, EINF , EO , Z , ETA ) ! DEFINES SET OF EQS 

DOUBLE PRECISION EK ( 0 : 4 , 0 : 1 0 ) , YA ( 0 : 1 0 ) , H , XA , PARTl ( 0 : 4 , 0 : 1 0 ) 

DOUBLE PRECISION PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTl ( J, 1 )=DEXP( (2.*Z**2.)/( (XA+1 . ) **2 . ) )/( ( XA+1 . ) **3 . ) 

PART2 ( J,1 )=( DEXP( ( Z**2 . )/( (XA+1 . ) **2 . ) ) ) * ( ETA*EINF/E0 ) 

EK ( J , 1 ) = ( 2 . *ETA-PART2 (J,l)*YA(l)-( ETA* ( EINF/EO ) * *2 . ) *PARTl ( J , 1 ) ) *H 
EK(J,2)= the second ode, if nec. 

RETURN 

END 
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n n 



THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



DOUBLE PRECISION YA ( 0 : 1 0 ) , YN ( 0 ; 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 
DOUBLE PRECISION XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 
DOUBLE PRECISION NI , NE , NI 1 , NI 2 , NI 3 , NI 2A , NI 3A , NEl 
REAL PI,XL,XCHK,EINF,EO,Z,ETA 
REAL EPSO,CHRG,K,TEMP,NINF,A 
INTEGER NS, PI 

EPS0=8 . 853742E-12 
CHRG=1 . 602E-19 
K=1.38054E-23 

PRINT * 

PRINT FOURTH ORDER RUNGE-KUTTA SCHEME ' 

PRINT FOR THE CARTESIAN COMPUTATION ' 

PRINT OF NI AND NE ' 

PRINT *, 'INPUT VALUES FOR EINF, EO , Z, AND ETA' 

READ *, EINF, E0,Z, ETA 

PRINT *, 'INPUT VALUES FOR TEMP. AND NINF' 

READ *, TEMP, NINF 

PRINT *, 'INPUT VALUE FOR a' 

READ *,A 

1 IM=1 ! Number of equations 

Y(l) = 1.00 ! Initial condition for y~l at x=XP. 

C Y(2) = 0 ! Initial condition for y~2 at x=XP (if nec) 

PRINT * 

PRINT *, 'INTERVAL OF X FOR PRINTING ?' 

READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-x)' 

READ *,H 

NS = NINT(PI/H) 

PRINT *, 'MAXIMUM X TO STOP CALCULATION ?' 

READ *,XL 

PRINT * , ' H= ' ,H 

P1 = 0 
XP=0 
HH=H/2 
KPLUS=1 . 0 



print * , ' HH= ' ,HH 
print * , ' NS= ' , NS 
PRINT * 



LI = 0 ! Line no. initialization 

PRINT*, 'LINE y/a KPLUS' 

WRITE (*,98) LI,XP, KPLUS 



74 



LI=LI+1 
DO N=1,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 



! Old time 
! New time 
! Midpoint time 



J=1 ! This part computes k~l. 

DO 1=1, IM 

YA( I )=Y( I ) 

END DO 
XA=XB 

CALL FUNCT(EK, J, YA,H,XA,EINF,EO,Z,ETA) 

J=2 ! This part computes k~2. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END DO 
XA=XM 

CALL FUNCT ( EK , J , YA , H , XA , El NF , EO , Z , ETA ) 

J=3 ! This part computes k~3. 

DO 1=1, IM 

YA(I)=Y(I)+EK(2,I)/2 
END DO 
XA=XM 

CALL FUNCT(EK, J, YA,H,XA,EINF,EO,Z,ETA) 

J=4 I This part computes k~4. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 3, I ) 

END DO 
XA=XP 

CALL FUNCT(EK, J,YA,H,XA,EINF,EO,Z,ETA) 



DO 1=1, IM I 4-th order Runge-Kutta scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2 , I ) * 2+EK ( 3 , I ) * 2+EK ( 4 , I ) ) /6 
END DO 



W=Y( 1 ) 

You now have W( y/a ) — to get K+( y/a ) you must 
multiply by { E( y/a ) / E(0) } 

E=(EINF/EO)*DEXP( (Z**2. )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

To get the derivative of W [ W' ( y/a ) ], just plug the 
computed values of W back in the original Ist-order ODE 



J=5 

YA ( 1 ) =W 
XA=XP 

CALL FUNCT(EK, J,YA,H,XA,EINF,EO,Z,ETA) 

WPRM=EK( 5,1 )/H 

NOW, to get d/dy [ K+ ] from d/dy [ K+/(E/E0) ] , must 

perform the following operation: 

KPLSPRM=(WPRM* ( E**2 . ) +KPLUS*E* ( -2 . * ( Z* * 2 . )/( (XP+1. )**3. ) ) )/E 



75 



c 

c 

c 



c 

c 

c 

c 

72 

73 

74 

75 

44 

88 

98 

200 



AND NOW, to generate the non-dimensionalized n-curves 
ne( y/a ) and ni( y/a ) 
just compute these equations; 

NIl=( KPLUS/E ) * ( EINF/EO ) 

Nl2A=(EPS0*K*TEMP)/( (CHRG**2, ) *NINF* ( A**2 . ) ) 

NI2=NI2A*( (6.*(Z**2.)/( (XP+1. )**4.))+(4.*(Z**4.)/( (XP+1. )**6. ) ) ) 

NI3A=( ( EPS0*E0 )/( CHRG*NINF*A) ) *E 
NI3=NI3A*(-2.*(Z**2. )/( (XP+1. )**3. ) ) 

NI=. 5*(NI1+NI2+NI3 ) 

NE1=( ( E0*EPS0* ( Z**2 . ) )/( CHRG*NINF*A) ) *E* ( 2 ./( (XP+1 . ) **3 . ) ) 
NE=NI+NE1 



To prevent every data point from being written to the file 
(resulting in unplottable 50,000 pt. data files), the following 
edits out a percentage of the data depending on 'where' it 
was generated. 



XCHK=(XP/H) 



,GT, 

,GT, 

,GT, 

,GT, 



100000) GOTO 72 
10000) GOTO 73 
1000) GOTO 74 
100) GOTO 75 



IF (XCHK 
IF (XCHK 
IF (XCHK 
IF (XCHK 
GOTO 44 
P1=P1+1 

IF (PI .NE. 10000) GOTO 88 
P1 = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi 
Pl = 0 



,NE. 1000) GOTO 88 



,NE. 100) GOTO 88 



.NE. 10) GOTO 88 



WRITE (10,*) XP, NI 
WRITE (11,*) XP, NE 

CONTINUE 



END DO 

WRITE (*,98) LI,XP, KPLUS 
FORMAT(lX, 12, FlO.6, 2X, 1P4E16.8) 

IF (XP .LT. XL) GOTO 28 

PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ * , Kl 

IF(Kl.EQ.l) GOTO 1 
PRINT* 

76 



END 



'k'k'k'k'f^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k 



SUBROUTINE FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA ) I DEFINES SET OF EQS 
DOUBLE PRECISION EK ( 0 : 4 , 0 : 10 ) , YA ( 0 : 1 0 ) , H , XA , PARTI ( 0 : 4 , 0 : 1 0 ) 

DOUBLE PRECISION PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTI ( J , 1 )=DEXP( (2.*Z**2.)/( (XA+1 . ) **2 . ) )/( (XA+1 . ) **3 . ) 

PART2 ( J , 1 )=( DEXP( ( Z**2 . )/( (XA+1 . ) **2 . ) ) ) * ( ETA*EINF/E0 ) 

EK( J, 1 )=( 2 . *ETA-PART2 ( J, 1 ) *YA( 1 )-( ETA* ( E INF/E 0 ) **2 . ) * PARTI ( J, 1 ) ) *H 
EK(J,2)= the second ode, if nec. 

RETURN 

END 
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n n 



THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



1 

C 



DOUBLE PRECISION YA( 0 : 10 ) , YN( 0 : 10 ) , EK( 0 : 4 , 0 : 10 ) , Y( 0 : 10 ) , XA 
DOUBLE PRECISION XP , XB , XM , H , HH , E , W, KPLUS , WPRM , KPLSPRM 
REAL PI,XL,XCHK,EINF,ERO,Q,ETA 
INTEGER NS, Pi 
PRINT * 

FOURTH ORDER RUNGE-KUTTA SCHEME ' 

FOR THE CYLINDRICAL COMPUTATION ' 

OF W, KPLUS, & KPLSPRM' 

( NORMAL-SIZED ANODE: rO = 1 cm ) ' 



PRINT 

PRINT 

PRINT 

PRINT 



* • 

t 

* > 

t 

■k • 



PRINT *, 'INPUT VALUES FOR EINF, ErO, Q, AND ETA' 

READ *, EINF, ERO, Q, ETA 

IM=1 ! Number of equations 

Y(l) = 1.00 I Initial condition for wl at r~=XP. 

Y(2) = 0 ! Initial condition for w2 at r~=XP (if nec) 



PRINT * 

PRINT *, 'INTERVAL OF r~ FOR PRINTING ?' 
READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-r~)' 
READ *,H 

NS = NINT(PI/H) 

PRINT *, 'MAXIMUM r~ TO STOP CALCULATION ?' 
READ *,XL 

PRINT * , ' dr~ = ' , H 

Pl = 0 
XP=0 
HH=H/2 . 

KPLUS=1 . 0 



print * , ' NS= ' , NS 
PRINT * 



28 



C 



LI = 0 

PRINT*, 'LINE r~ 
WRITE (*,98) LI,XP, 



! Line no. initialization 
KPLUS' 

KPLUS 



LI=LI+1 
DO N=1,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 



! Old time 
! New time 
! Midpoint time 



Runge-Kutta Scheme 



J = 1 

DO 1=1, IM 



This is part 1. 
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YA( I )=Y( I ) 

END DO 
XA=XB 

CALL FUNCT(EK,J 
J = 2 

DO 1=1, IM 

YA( I )=Y( I )+EK 
END DO 
XA=XM 

CALL FUNCT(EK,J 
J = 3 

DO 1=1, IM 

YA( I )=Y( I )+EK 
END DO 
XA=XM 

CALL FUNCT(EK 
J = 4 

DO 1=1, IM 

YA( I )=Y( I )+EK 
END DO 
XA=XP 

CALL FUNCT(EK,J 

DO 1=1, IM 

Y(I)=Y(I)+(E 
END DO 



YA,H,XA,EINF,ERO,Q 
I This is part 2. 

1.1) /2 

YA,H,XA,EINF,ERO ,Q 
! This is part 3. 

2.1) /2 

J,YA,H,XA,EINF,ERO 
! This is part 4. 

3.1) 

YA,H,XA,EINF,ERO,Q 

! 4-th order Runge 
(1,1 )+EK( 2 , I ) *2+EK 



ETA) 



ETA) 



Q, ETA) 



ETA) 

Kutta scheme 
3 , 1 ) *2+EK( 4 , 1 ) )/6 



W=Y( 1 ) 

You now have W( r~ ) — to get K+( r~ ) you must 
multiply by { E( r~ ) / E(0) } 

E=( EINF/ERO ) *DEXP( (Q**2. )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

To get the derivative of W [ W' ( r~ ) ] < just plug the 

computed values of W back in the original Ist-order ODE 



J=5 

YA( 1 )=W 
XA=XP 

CALL FUNCT ( EK , J , YA , H , XA, EINF , ERO , Q , ETA ) 

WPRM=EK( 5, 1 )/H 

NOW, to get d/dr~ [K+] from d/dr~ [ K+/(E/ErO) ] , must 

perform the following operation: 

KPLSPRM=(WPRM* ( E**2 . )+KPLUS*E* ( -2.*(Q**2.)/( (XP+1 . ) **3 . ) ) )/E 

To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on 'where' they occur 

XCHK=(XP/H) 

IF (XCHK .GT. 100000) GOTO 72 jq 

IF (XCHK .GT. 10000) GOTO 73 



IF (XCHK .GT. 1000) GOTO 74 
IF (XCHK .GT. 100) GOTO 75 
GOTO 44 



72 

73 

74 

75 



Pl=Pl+l 
IF (PI .NE. 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi .NE. 
P1 = 0 
GOTO 44 
P1=P1+1 
IF (Pi .NE. 
P1 = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi .NE. 
Pl = 0 



10000) GOTO 88 



1000) GOTO 88 



100) GOTO 88 



10) GOTO 88 



44 


WRITE 


(13,*) 


XP, 


KPLSPRM 




WRITE 


(12,*) 


XP, 


KPLUS 




WRITE 


(14,*) 


XP, 


W 



88 CONTINUE 



END DO 

WRITE (*,98) LI,XP, KPLUS 
98 FORMATdX, I2, FlO.6, 2X, 1P4E16.8) 

IF (XP .LT. XL) GOTO 28 

200 PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ *,K1 

IF(Kl.EQ.l) GOTO 1 
PRINT* 

END 



SUBROUTINE FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) ! DEFINES SET OF EQS 

DOUBLE PRECISION EK ( 0 : 4 , 0 : 1 0 ) , YA( 0 : 1 0 ) , H , XA , PARTl ( 0 : 4 , 0 : 10 ) 

DOUBLE PRECISION PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTI ( J,1 )=ETA*(DEXP( (2.*Q**2. )/( (XA+1. )**2. ) )/( (XA+1. )**3. ) ) 
PART2( J,1)=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) ) * ( ETA*EINF/ER0 ) 

EK( J,1)=(2.*ETA-PART2( J,l)*YA(l)-( ( EINF/ERO ) * *2 . )*PARTl( J,l) ) *H 
C EK(J,2)= the second ode, if nec. 

RETURN 

END 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



DOUBLE PRECISION YA( 0 : 1 0 ) , YN( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 

DOUBLE PRECISION XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 

DOUBLE PRECISION NI , NE , NI 1 , NI 2 , NI 3 , NEl 

REAL CCl ,CC2,CC3,S2,S3,S2A,S2B,S2C 

REAL PI ,XL,XCHK,EINF,ERO ,Q,ETA 

REAL EPSO,CHRG,K,TEMP,NINF,B,RO 

INTEGER NS, Pi 

EPS0=8 . 853742E-12 
CHRG=1 . 602E-19 
K=1 . 38054E-23 
R0= . 01 

PRINT * 

PRINT * , ' 

PRINT * , ' 

PRINT * , ' 

PRINT * , ' 



FOURTH ORDER RUNGE-KUTTA SCHEME ' 
FOR THE CYLINDRICAL COMPUTATION ' 

OF NE AND NI' 

( NORMAL-SIZED ANODE: rO = 1 cm ) ' 



PRINT *, 'INPUT VALUES FOR EINF, ErO, Q, AND ETA' 

READ *, EINF, ERO, Q, ETA 

PRINT *, 'INPUT VALUES FOR TEMP. AND NINF' 

READ *, TEMP, NINF 

PRINT *, 'INPUT VALUE FOR b' 

READ * , B 

IM=1 ! Number of equations 

Y(l) = 1.00 I Initial condition for wl at r~=XP. 

Y(2) = 0 ! Initial condition for w2 at r~=XP (if nec) 

PRINT * 

PRINT *, 'INTERVAL OF r~ FOR PRINTING ?' 

READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-r~)' 

READ *,H 

NS = NINT(PI/H) 

PRINT *, 'MAXIMUM r" TO STOP CALCULATION ?' 

READ *,XL 

PRINT *, ' dr" = ' ,H 

Pl = 0 
XP = 0 
HH=H/2 . 

KPLUS=1 . 

print * , ' NS= ' , NS 
PRINT * 



LI = 0 ! Line no. initialization 

PRINT*, 'LINE r~ KPLUS' 

WRITE (*,98) LI,XP, KPLUS 
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28 



C 



C 



LI=LI+1 
DO N=1,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 



! Old time 
! New time 
! Midpoint time 



Runge-Kutta Scheme 



J=1 ! This is part 1. 

DO 1=1, IM 

YA(I)=Y(I) 

END DO 
XA=XB 

CALL FUNCT( EK, J, YA,H,XA,EINF,ERO,Q,ETA) 

J=2 I This is part 2. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END DO 
XA=XM 

CALL FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) 

J=3 I This is part 3. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 2 , I )/2 
END DO 
XA=XM 

CALL FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) 

J=4 ! This is part 4. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 3 , I) 

END DO 
XA=XP 

CALL FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) 

DO 1=1, IM ! 4-th order Runge-Kutta scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2 , I ) *2+EK( 3 , I ) *2+EK( 4 , I ) )/6 
END DO 



W=Y( 1 ) 

c You now have W( r~ ) — to get K+( r~ ) you must 

c multiply by { E( r~ ) / E(0) } 

E=( EINF/ERO ) *DEXP( (Q**2. )/( (XP+1 . )**2. ) ) 

KPLUS=W*E 

c To get the derivative of W [ W' ( r~ ) ] , just plug the 

c computed values of W back in the original Ist-order ODE 



J=5 

YA( 1)=W 
XA=XP 

CALL FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) 
WPRM=EK( 5,1)/H 09 



NOW, to get d/dr~ [ K+ ] from d/dr~ [ K+/(E/ErO) ] , must 

perform the following operation: 

KPLSPRM=(WPRM* ( E**2 . ) +KPLUS*E* ( -2 . * ( Q**2 . )/( (XP+1. )**3. ) ) )/E 



AND NOW, to generate the nondimensionalized n-curves 
ne( r~ ) and ni( r~ ) 
compute the following equations 

NI1=( KPLUS/E ) * ( EINF/ERO ) 

CC2=( K*TEMP*EPSO )/( ( CHRG**2 . ) *NINF* ( B**2 . ) ) 

S2A=( 6 . * ( Q**2 . ) )/( (XP+1. )**4. ) 

S2B=(4.*(Q**4.))/( (XP+1. )**6. ) 

S2C=( 1/(XP+(R0/B) ))*((-2.*(Q**2.))/( (XP+1. )**3. ) ) 
S2=S2A+S2B+S2C 
NI2=CC2*S2 

CC3=( EPS0*ER0 )/( CHRG*NINF*B) 

S3=( (-2.*(Q**2. ) )/( (XP+1. )**3. ) )+(l./(XP+(R0/B) ) ) 
NI3=CC3*E*S3 

NI=.5*(NI1+NI2+NI3) 

CC1=( EPS0*ER0 )/( CHRG*NINF*B ) 

NE1=CC1*E*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) + ( 1 . /( XP+ ( RO/B ) ) ) ) 
NE=NI-NEl 



To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on 'where' they occur 



XCHK=(XP/H) 

IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
GOTO 44 



100000) GOTO 72 
10000) GOTO 73 
1000) GOTO 74 
100) GOTO 75 



Pl=Pl+l 
IF (Pi .NE. 
Pl = 0 
GOTO 44 
P1=P1+1 
IF (Pi .NE. 
Pl = 0 
GOTO 44 
P1=P1+1 
IF (Pi .NE. 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (PI .NE. 
Pl = 0 



10000) GOTO 88 



1000) GOTO 88 



100) GOTO 88 



10) GOTO 88 



WRITE (10,*) XP, NI 
WRITE (11,*) XP, NE 



CONTINUE 
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END DO 



WRITE (*,98) LI,XP, KPLUS 
98 FORMATdX, I2, FlO.6, 2X, 1P4E16.8) 

IF (XP .LT. XL) GOTO 28 

200 PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ *,K1 

IF(Kl.EQ.l) GOTO 1 
PRINT* 

END 

Qicic’k'k'k'k'k'k'k'k’k’k’k'k'k’k’k’k'k'kic'k'k'k'k'k'kic'k'k'k'k'k'k'k'k'k'k'k'k'k 



SUBROUTINE FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) ! DEFINES SET OF EQS 

DOUBLE PRECISION EK ( 0 : 4 , 0 : 10 ) , YA( 0 : 1 0 ) , H , XA , PARTl ( 0 : 4 , 0 : 1 0 ) 

DOUBLE PRECISION PART2 ( 0 : 4 , 0 : 1 0 ) 

PART1( J, 1 )=ETA*( DEXP( ( 2 . *Q**2 . )/( (XA+1 . ) **2 . ) )/( (XA+1 . ) **3 . ) ) 
PART2( J,1) = (DEXP( (Q**2. )/( (XA+1. )**2. ) ) ) * ( ETA*EINF/ER0 ) 

EK( J,1)=(2.*ETA-PART2( J,l)*YA(l)-( ( EINF/ERO ) * * 2 . )*PARTl( J,l) ) *H 
C EKd,2)= the second ode, if nec. 

RETURN 

END 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



DOUBLE PRECISION YA ( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 ; 1 0 ) , Y ( 0 ; 1 0 ) , XA 
DOUBLE PRECISION XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 
REAL PI ,XL,XCHK,EINF,ER0,Q,ETA,R0,G,B,EPS0,CHRG 
INTEGER NS, Pi 



PRINT 


* 










PRINT 




FOURTH 


ORDER RUNGE- 


KUTTA 


SCHEME ' 


PRINT 




FOR THE CYLINDRICAL 


COMPUTATION ' 


PRINT 


★ f 

f 


OF 


W, KPLUS, & 


KPLSPRM' 


PRINT 


-k f 

f 


( WIRE- 


THIN ANODE: 


rO = 


0.1 mm ) ' 



PRINT 
READ * 


* , ' INPUT VALUES 
, EINF, ERO ,Q, ETA 


FOR 


EINF, 


ErO, Q, AND ETA' 


PRINT 
READ * 


* , ' INPUT VALUES 
,G,B 


FOR 


GAMMA 


AND B' 



R0=.0001 

EPS0=8 .853742E-12 
CHRG=1 .602E-19 

IM=1 1 Number of equations 

Y(l) = 1.00 ! Initial condition for wl at r~=XP. 

Y(2) = 0 ! Initial condition for w2 at r~=XP (if nec) 

PRINT * 

PRINT *, 'INTERVAL OF r~ FOR PRINTING ?' 

READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-r~)' 

READ *,H 

NS = NINT(PI/H) 

PRINT *, 'MAXIMUM r~ TO STOP CALCULATION ?' 

READ *,XL 

PRINT * , ' dr~ = ' , H 

Pl = 0 
XP=0 
HH=H/2 . 

KPLUS=1 . 

print * , ' NS= ' , NS 
PRINT * 



LI = 0 I Line no. initialization 

PRINT*, 'LINE r~ KPLUS' 

WRITE (*,98) LI,XP, KPLUS 



LI=LI+1 
DO N=1,NS 
XB=XP 
XP=XP+H 



! Old time 
! New time 



n n 



XM=XB+HH ! Midpoint time 



C Runge-Kutta Scheme 

J=1 ! This is part 1. 



DO 1 = 1, IM 

YA( I )=Y( I ) 

END DO 
XA=XB 

CALL FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) 

J=2 ! This is part 2. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END DO 
XA=XM 

CALL FUNCT(EK, J,YA,H,XA, EINF, ERO, Q, ETA, G,B) 

J=3 ! This is part 3. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 2, I )/2 
END DO 
XA=XM 

CALL FUNCT(EK, J,YA,H,XA, EINF, ERO, Q, ETA, G,B) 

J=4 ! This is part 4. 

DO 1 = 1, IM 

YA( I )=Y( I )+EK( 3 , I ) 

END DO 
XA=XP 

CALL FUNCT(EK, J,YA,H,XA, EINF, ERO, Q, ETA, G,B) 

DO 1=1, IM I 4-th order Runge-Kutta scheme 

Y( I )=Y( I)+( EK( 1 , I )+EK( 2 , I ) *2+EK( 3 , I ) *2+EK( 4, I ) )/6 
END DO 



W=Y( 1 ) 

c You now have W( r~ ) — to get K+( r~ ) you must 

c multiply by { E( r~ ) / E(0) } 

E=( EINF/ERO ) *DEXP( (Q**2. )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

c To get the derivative of W [ W' ( r~ ) ] , just plug the 

c computed values of W back in the original Ist-order ODE 



J=5 

YA( 1 )=W 
XA=XP 

CALL FUNCT ( EK,J ,YA, H,XA, EINF, ERO ,Q, ETA, G,B ) 

WPRM=EK( 5 , 1 )/H 

NOW, to get d/dr~ [ K+ ] from d/dr~ [ K+/(E/ErO) ] , must 

perform the following operation: 

KPLSPRM=( WPRM* ( E**2 . ) +KPLUS*E* ( -2 . * ( Q**2 . )/( (XP+1. )**3. ) ) )/E 
C To keep from generating unplottable 50,000 point data files. 
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the following edits out some of the data points depending 
on 'where' they occur 



XCHK=(XP/H) 

IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
GOTO 44 



100000) GOTO 72 
10000) GOTO 73 
1000) GOTO 74 
100) GOTO 75 



P1=P1+1 
IF (Pi .NE. 
P1 = 0 
GOTO 44 
Pl=Pl+l 
IF (PI .NE. 
Pl = 0 
GOTO 44 
Pl=Pl+l 



10000) GOTO 88 



1000) GOTO 88 



IF (Pi .NE. 100) GOTO 88 
P1 = 0 
GOTO 44 
P1=P1+1 



IF (PI 


.NE. 


10) GOTO 88 


Pl = 0 










WRITE 


(13, 


*) 


XP, 


KPLSPRM 


WRITE 


(12, 


*) 


XP, 


KPLUS 


WRITE 


(14, 


*) 


XP, 


W 



CONTINUE 



END DO 

WRITE (*,98) LI,XP, KPLUS 
FORMATdX, 12, F10.6, 2X, 1P4E16.8) 

IF (XP .LT. XL) GOTO 28 

PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ *,K1 

IF(Kl.EQ.l) GOTO 1 
PRINT* 

END 



SUBROUTINE FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) ! DEFINES ODEs 

DOUBLE PRECISION EK ( 0 : 4 , 0 : 1 0 ) , YA ( 0 : 1 0 ) , PART2 ( 0 : 4 , 0 : 1 0 ) 

DOUBLE PRECISION PARTS ( 0 : 4 , 0 : 1 0 ) , PART4 ( 0 : 4 , 0 : 1 0 ) , H , XA , Cl 

REAL RO 

R0=.0001 

C2=ETA*G* ( ( EINF/ERO ) **2 . ) 

PART4( J,1 )=C2*(DEXP( ( 2 . * ( Q* * 2 . ) ) /( (XA+1. )**2. ) ) )/(XA*B+R0) 
C1=ETA*( (G/R0)+1. )*( ( EINF/ERO ) **2 . ) 

PART3( J,1 )=C1*(DEXP( (2.*Q**2. )/( (XA+1. )**2. ) )/( (XA+1. )**3. ) ) 
PART2( J,1 )=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) )*( ETA* EINF/ERO ) 

EK ( J , 1 ) = ( 2 . *ETA-PART2 ( J , 1 ) * YA ( 1 ) -PARTS ( J , 1 ) +PART4 ( J , 1 ) ) *H 
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C EK(J,2)= the second ode, if nec. 

RETURN 
END 
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THIS PROGRAM IS BUILT AROUND A GRAFkit 3.1 RUNGE-KUTTA 
ALGORITHM PRODUCED BY SCO, INC. 



DOUBLE PRECISION YA( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 

DOUBLE PRECISION XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 

DOUBLE PRECISION NI , NE , NI 1 , NI 2 , NI 3 , NEl 

REAL CCl , CC2 ,CC3,S2,S3,S2A,S2B,S2C,S2D 

REAL PI ,XL,XCHK,EINF,ERO,Q,ETA 

REAL EPSO , CHRG , K , TEMP , NI NF , B , G , RO 

INTEGER NS, Pi 

EPS0=8 . 853742E-12 
CHRG=1.602E-19 
K=1.38054E-23 
R0=.0001 

PRINT * 

PRINT * , ' 

PRINT * , ' 

PRINT * , ' 

PRINT * , ' 



FOURTH ORDER RUNGE-KUTTA SCHEME ' 
FOR THE CYLINDRICAL COMPUTATION ' 
OF NE AND NI' 

( WIRE-THIN ANODE: rO = 0.1 mm )' 



PRINT *, 'INPUT VALUES FOR EINF, ErO, Q, AND ETA' 

READ *, EINF, ERO, Q, ETA 

PRINT *, 'INPUT VALUES FOR TEMP. AND NINF' 

READ *, TEMP, NINF 

PRINT *, 'INPUT VALUES FOR GAMMA AND b' 

READ * , G , B 

IM=1 ! Number of equations 

Y(l) = 1.00 ! Initial condition for wl at r~=XP. 

Y(2) = 0 ! Initial condition for w2 at r~=XP (if nec) 

PRINT * 

PRINT *, 'INTERVAL OF r~ FOR PRINTING ?' 

READ *,PI 

PRINT *, 'INPUT THE STEP SIZE (delta-r~)' 

READ *,H 



NS = NINT(PI/H) 

PRINT *, 'MAXIMUM r~ TO STOP CALCULATION ?' 
READ *,XL 

PRINT * , ' dr~ = ' ,H 



Pl = 0 
XP=0 
HH=H/2 . 
KPLUS=1 . 



print * , ' NS= ' , NS 
PRINT * 



LI = 0 ! Line no. initialization 

PRINT*, 'LINE r~ KPLUS' 

WRITE (*,98) LI,XP, KPLUS 
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28 LI=LI+1 

DO N=1,NS 

XB=XP ! Old time 

XP=XP+H ! New time 

XM=XB+HH ! Midpoint time 

C Runge-Kutta Scheme 

J=1 ! This is part 1. 

DO 1 = 1, IM 

YA(I)=Y(I) 

END DO 
XA=XB 

CALL FUNCT( EK, J, YA,H,XA,EINF,ER0,Q,ETA,G,B) 

J=2 ! This is part 2. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END DO 
XA=XM 

CALL FUNCT(EK, J,YA,H,XA,EINF,ERO,Q,ETA,G,B) 

J=3 ! This is part 3. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 2 , I )/2 
END DO 
XA=XM 

CALL FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) 

J=4 ! This is part 4. 

DO 1=1, IM 

YA( I )=Y( I )+EK( 3, I ) 

END DO 
XA=XP 

CALL FUNCT (EK, J,YA,H,XA, EINF, ERO, Q, ETA, G,B) 

DO 1=1, IM ! 4-th order Runge-Kutta scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2 , I ) *2+EK( 3,1) *2+EK( 4 , I ) ) /6 
END DO 

C 



W=Y( 1 ) 

c You now have W( r~ ) — to get K+( r~ ) you must 

c multiply by { E( r~ ) / E(0) } 

E=( EINF/ERO ) *DEXP( (Q**2. )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

c To get the derivative of W [ W'( r"' ) ], just plug the 

c computed values of W back in the original Ist-order ODE 



J=5 

YA( 1 )=W 
XA=XP 

CALL FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) 
WPRM=EK( 5 , 1 )/H 
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NOW, to get d/dr~ [ K+ ] from d/dr~ [ K+/(E/ErO) ] , must 

perform the following operation: 

KPLSPRM=(WPRM* ( E**2 . ) +KPLUS*E* ( -2 . * ( Q* * 2 . )/( (XP + 1. )**3. ) ) )/E 



AND NOW, to generate the nondimensionalized n-curves 
ne( r~ ) and ni( r~ ) 
compute the following equations 

NIl= ( KPLUS/E ) * ( EINF/ERO ) 

CC2=( K*TEMP*EPSO )/( ( CHRG**2 . ) *NINF* ( B**2 . ) ) 

S2A=( 6 . * ( Q**2 . ) )/( (XP+1 . ) **4 . ) 

S2B=( 4 . * (Q**4 . ) )/( (XP + 1 . ) **6 . ) 

S2C=(1./(XP+(R0/B) ) )*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) ) 
S2D=(l./( (XP+(R0/B) )**2. ) ) 

S2=S2A+S2B+S2C-S2D 

NI2=CC2*S2 

CC3=( EPS0*ER0 )/( CHRG*NINF*B) 

S3=( (-2.*(Q**2. ) )/( (XP+1. )**3. ) )+( 1 ./(XP+(R0/B) ) ) 
NI3=CC3*E*S3 

NI=.5*(NI1+NI2+NI3 ) 

CC1=( EPS0*ER0 )/( CHRG*NINF*B) 

NE1=CC1*E*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) + ( 1 . /( XP+ ( RO/B ) ) ) ) 
NE=NI-NE1 



To keep from generating unplottable 50,000 point data files, 
the following edits out some of the data points depending 
on 'where' they occur 



XCHK=(XP/H) 

IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
IF (XCHK .GT. 
GOTO 44 



100000) GOTO 72 
10000) GOTO 73 
1000) GOTO 74 
100) GOTO 75 



Pl=Pl+l 
IF (PI .NE. 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (PI .NE. 
P1 = 0 
GOTO 44 
P1 = P1 + 1 
IF (Pi .NE. 
Pl = 0 
GOTO 44 
Pl=Pl+l 
IF (Pi .NE. 
Pl = 0 



10000) GOTO 88 



1000) GOTO 88 



100) GOTO 88 



10) GOTO 88 



WRITE (10,*) XP, NI 
WRITE (11,*) XP, NE 
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88 



CONTINUE 



END DO 

WRITE (*,98) LI,XP, KPLUS 
98 FORMATdX, 12, FlO.6, 2X, 1P4E16.8) 

IF (XP .LT. XL) GOTO 28 

200 PRINT* 

PRINT*, 'TYPE 1 TO CONTINUE, OR 0 TO STOP.' 
READ *,K1 

IF(Kl.EQ.l) GOTO 1 
PRINT* 

END 



SUBROUTINE FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) I DEFINES ODES 
DOUBLE PRECISION EK ( 0 : 4 , 0 : 1 0 ) , YA( 0 : 1 0 ) , H , XA , PART2 ( 0 : 4 , 0 : 1 0 ) 
DOUBLE PRECISION PARTS ( 0 : 4 , 0 : 1 0 ) , PART4 ( 0 : 4 , 0 : 1 0 ) 

REAL R0,KK4,KK3 
R0=.0001 

KK4=ETA*G* ( ( EINF/ERO ) **2 . ) 

PART4 ( J, 1 )=KK4* ( DEXP ( (2.*(Q**2.))/( (XA+1 . ) **2 . ) ) )/(XA*B+R0 ) 
KK3=ETA*( (G/R0)+1. )*( ( EINF/ERO )** 2 . ) 

PART3( J,1)=KK3*(DEXP( (2.*Q**2. )/( (XA+1. )**2. ) )/( (XA+1. )**3. ) ) 
PART2( J,1) = (DEXP( (Q**2. )/( (XA+1. )**2. ) ) )*( ETA* EINF/ERO ) 

EK( J , 1 )=( 2 . *ETA-PART2( J, 1 ) *YA( 1 ) -PARTS ( J , 1 ) +PART4 ( J, 1 ) ) *H 
C EK(J,2)= the second ode, if nec. 

RETURN 

END 
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